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Abstract. Regge calculus is used to construct initial data for vacuum 
axisymmetric Brill waves at a moment of time symmetry. We argue that only 
a tetrahedral lattice can successfully reproduce the continuum solution, and 
develop a simplicial axisymmetric lattice based on the co-ordinate structure of 
the continuum metric. This is used to construct initial data for Brill waves in an 
otherwise flat spacetime, and for the distorted black hole spacetime of Bernstein. 
These initial data sets are shown to be second order accurate approximations to 
the corresponding continuum solutions. 



PACS numbers: 

1. Introduction 

Since it's inception in 1961, Regge calculus |Q has been extensively studied in highly 
symmetric spacetimes, for which corresponding exact solutions of Einstein's equations 
are often available. A complete review and bibliography of this early work is provided 
by Williams and Tuckey Q . 

Following the realization that a fully decoupled, parallelizable, three-plus-one 
dimensional evolution scheme occurs naturally within Regge calculus || ||, this 
simplicial approach to gravity is now on the verge of tackling physically interesting 
and dynamic problems. The first tentative step along this path was recently completed, 
with the successful application of simplicial Regge calculus to the Kasner spacetime 
in (3 + l)-dimensions ||. 

A vital precursor to the evolution problem in any general relativistic simulation 
is the construction of consistent initial data. Gentle and Miller j|] present a general 
prescription for the calculation of two-surface initial data for the Regge lattice, 
although their approach has only been applied to the Kasner cosmology. As a prelude 
to the construction of simplicial initial data for complex spacetimes, in this paper we 
consider the restricted case of vacuum, axisymmetric initial data at a moment of time 
symmetry. This is the first step towards our goal, and provides a benchmark against 
which future, fully four-dimensional, initial data may be compared. 

We construct vacuum initial data for the pure Brill wave spacetime j7|, and for 
Brill waves in a black hole spacetime — the "distorted black holes" first considered 
by Bernstein Standard finite-difference techniques have been used to study pure 
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Brill wave spacetimes by Eppley ||, Miyama (T^|, Holz et al |TlJ and Alcubierre et 
al | p^ |. The axisymmetric Brill wave plus black hole spacetime has been extensively 
investigated by Bernstein and more recently, a generalization to (3 + l)-dimensions 
has been considered by Camarda jl3| . 

Dubai Q constructed pure Brill wave initial data using Regge calculus, with 
a lattice based on prisms rather than simplices. A similar approach was taken 
in an earlier version of the current work [ fi"5| . We will consider the advantages 
and disadvantages of using a prism-based lattice, and argue that the prism-based 
construction cannot, in general, satisfactorily reproduce the continuum solution, even 
in axisymmetry. 

We begin in section || with a short survey of the important results from the 
continuum theory. In section [| we discuss the initial value problem at a moment 
of time symmetry in Regge calculus, and then proceed in section |i] to examine an 
axisymmetric lattice built from prisms. In section ^| we turn to a fully simplicial 
lattice, and finally, in section ^, we analyze the convergence of the simplicial initial 
data to the continuum solution. 



2. The time symmetric initial value problem 

The vacuum constraint equations of General Relativity may be written in the form 
(see, for example, Misner, Thorne and Wheeler |l6||) 

R+ {trK) 2 +K ab K ab = (1) 
V b (K ab - 7 afc trK) =0 (2) 

where R a bcd is the intrinsic curvature tensor on the spacelike hypersurface, and R is 
the Ricci scalar. The extrinsic curvature K ab determines the embedding of this slice 
in the spacetime. At a moment of time symmetry K ao = 0, and the constraints reduce 
to vanishing of the scalar three-curvature: 

R = 0. (3) 

The standard approach to the solution of this equation involves a conformal 
decomposition of the three-metric, such that 

lab = i> 4 Jab (4) 

which recasts equation (^) in the form 

VV = ~^R (5) 

which is linear in the conformal factor The base metric i ab is used to calculate the 
base scalar curvature R. The base three-geometry may be freely specified, although 
the choice may restrict the class of possible solutions. Solving this single linear 
elliptic partial differential equation yields the full initial data set at a moment of 
time symmetry. 



2.1. Pure gravitational radiation 

In this section we introduce the axisymmetric vacuum initial value problem posed 
by Brill [Q. Conformal decomposition using a flat base metric leads only to trivial 
solutions. The approach taken by Brill j?j was to introduce a metric of the form 

ds 2 = i> i {e 2q (dp 2 + dz 2 )+p 2 d0 2 }, (6) 




Figure 1. The conformal factor ip for Brill wave initial data, calculated on a 
601 X 601 grid, using a centred finite difference approximation to equation (R). 
The Eppley form of q(p, z) (equation (H)) was used, with wave amplitude a = 10, 
and the outer boundaries were placed at p = 20 and z = 20. 



where the arbitrary function q(p, z) can be considered the distribution of gravitational 
wave amplitude fT^] , Brill showed that q(p, z) must satisfy the boundary conditions 

q(0,z) = 0, q p (Q,z)=Q, and q z (p,0)=Q, 

together with the condition that q = 0(r~ 2 ) as r — » oo, to ensure that the hypersurface 
has an asymptotically well defined mass. With this choice of background metric, the 
Hamiltonian constraint (||) takes the form 

' d 2 q d 2 q s 



.t 

4 \dp 2 ' dz 



which is solved for ip(p, z) once q{p, z) is given. 

To allow comparison with Eppley ||, we choose q(p, z) to be of the form 

^2 



ap 



where 



P 



(7) 



(8) 



1 + r n 

and the boundary conditions on q imply that n > 4. In the remainder of this paper 
we set n = 5. The single remaining parameter, the wave amplitude a, is arbitrary. 

To obtain an axisymmetric, asymptotically flat solution which is reflection 
symmetric on the z = axis, we use boundary conditions on ip of the form 



dp 



§Vo) = o, 

oz 



(9) 



as 



(10) 



together with a Robin outer boundary condition, 

dip 1 — ip 
dr r 

Figure |l|a shows the solution of equation (0) obtained using a centred finite 
difference approximation to both the equation and boundary conditions. The solution 
was calculated on a 601 x 601 grid with a Brill wave amplitude a = 10. 
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2.2. Distorted black hole initial data 

We now turn to the Brill wave plus black hole spacetime, investigated in great detail 
by Bernstein ||. As before, a perturbation is introduced onto the background metric, 
and the conformal factor is calculated from the single initial value equation (||). Note 
that we solve the full initial value problem for the wave and black hole together. This 
is not a perturbation solution. 



Mirroring the pure Brill wave spacetime of section 2.1, we write the physical 
metric on the initial surface in the form 

ds 2 = ip A {e 2q (dn 2 + dO 2 ) + sin 2 9d<j) 2 } , (11) 

where the exponential radius co-ordinate r\ is defined from p — mexp(r])/2. The 
"mass" m is that of the black hole alone, obtained by setting q(n, 9) = 0. In (77, 9) 
co-ordinates, the isolated black hole solution takes the form 

ip bh = V?m cosh (77/2) . (12) 

The Hamiltonian constraint (||) evaluated on the black hole plus Brill wave metric, 
equation (^), yields the linear constraint 

2 , "0 ( d 2 q , d 2 q 



w—nw+w- 1 ) (i3) 

which is again solved for the conformal factor ip once the perturbation q(rj, 9) is given. 
We choose boundary conditions of the form 

g (M) _0, % M -0, (») 

on the "inner" boundaries, where rj = corresponds to an Einstein-Rosen bridge at 
r = m/2. These represent reflection symmetry on the 9 = 0, 9 = n/2 and r\ = 
boundaries. The outer boundary condition is again applied using the Robin approach, 
adapted to the exponential radial co-ordinate r\. Examining the asymptotic form of 
the black hole solution yields the condition 

**4.* = ./!2,1 - 



a, V2-"' < 15 » 

on the outer boundary. 

As a test of the axisymmetric code, we solve equation ([l3]) using a second-order 
finite difference approximation with 3(77, 9) — 0, together with the boundary conditions 
above. We expect that the solution tpij will converge as the second power of the 
grid spacing to the analytic black hole solution, given by equation (111). We have 
performed such a calculation, and find that the axisymmetric numerical solution does 
indeed converge as the second power of the grid spacing to the black hole solution. 
This provides some confidence in the the continuum code, which is vital since it will 
be used later to benchmark the Regge calculations. 

Returning to the construction of distorted black hole initial data, we follow 
Bernstein and choose q(n, 9) to be of the form 

q = asm 2 9\e.- g + + e" s - | (16) 

with g± — (rj ± b) /u>, which contains the set of free parameters (a, b, to). 

The solution of the Hamiltonian constraint ( |l3| ) is shown in figure ^| for several 
illustrative choices of the free parameters. The figure shows the ratio of the conformal 
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(c) (d) 

Figure 2. The continuum solution for the black hole plus Brill wave spacetime. 
The conformal factor is plotted as a ratio of the black hole solution (ip/4>bh), with 
the results calculated on a 601 X 601 grid with m = 1. An artificial Cartesian 
projection is used, and this is defined in the text. The choice of parameters are 
(a) (1,0,1), (b) (-1,0,1), (c) (1,2,1) and (d) (-1,2,1). 

factor to the black hole solution, ip/ipbh- The outer boundary was placed at r\ = 6, 
and the results are plotted on a pseudo-Cartesian grid defined by 

x = (r\ + tjq) sin 9 

V = (»7 + ?7o)cos0 

where i]q = k is used to display an artificial throat in (77, 9) coordinates. These results 
compare well with the calculations of Bernstein ||, who also used a second order 
finite-difference approximation to the Hamiltonian constraint. 

3. The Regge initial value problem at a moment of time symmetry 

We found in section || that the continuum initial value problem at a moment of 
time symmetry reduces to the vanishing of the scalar three-curvature intrinsic to the 
hypersurface. By using this knowledge, the full two-slice Regge initial value problem 
can be reduced to a purely three-dimensional calculation. 

The Regge equivalent of the scalar curvature for a three dimensional hypersurface 
has been given by Wheeler |l7j . The Regge equivalent of R = is 

J2L ab Sab = (17) 
b 



where L ab is the edge joining vertex a to vertex b, and s a b is the deficit angle 
(curvature) about L ab . The summation is over all edges L a i, which meet at vertex 
a. 

In order to solve this restricted initial value problem, we mirror the continuum 
formulation and perform a conformal decomposition. This reduces the problem to the 
solution of the single Regge-Hamiltonian constraint, equation ([l?]), for the simplicial 
conformal factor. The conformal factor is defined on the vertices of the lattice, and we 
perform the decomposition on the lattice edge lengths, since they correspond closely to 
the continuum metric. For the edge L a b connecting vertices a and b, the decomposition 
takes the form 

Lab — tPab Lab 

where L a b is the freely chosen base edge length. Since the conformal factor is defined 
on the vertices, we use the second order accurate expression 

^ai = j W»a + Vb) (18) 

to define the conformal factor ip a b acting on the edge L a b, where ip a is the conformal 
factor at vertex a. 

So far we have considered only the restricted case of initial data at a moment of 
time symmetry — a purely three dimensional problem. However, it is possible to use 
the solution to equation ( p~7| ) to obtain full axisymmetric two-surface simplicial initial 
data. 

After solving equation (0), we have a set of physical edges {L a b}, which fix the 
geometry on a simplicial three-surface. Using the time-symmetry condition, we can 
construct two-slice initial data by creating two identical three-dimensional lattices 
using these edges, and then filling in the intervening four-volume with four-simplices. 
This involves introducing a timelike edge to join equivalent vertices in the two three- 
geometries, together with a "brace" edge which stretches from one surface to the next, 
with one such brace lying above each edge in the lower surface. This constructs a set 
of four-simplices between the two surfaces. Care must be taken in the choice of these 
braces if a particular lattice structure is desired. In particular, if we require a lattice 
amenable to the Sorkin evolution scheme ||. 

Assuming symmetry about the centre of the intervening volume, the squared 
lengths b 2 ab of the brace edges stretching from one surface to the next are 

b l+b = L lb + T a+a 

where vertex a + lies on the upper surface, and t 2 +q = r 2 < is the proper time 
between the two slices, measured along the timelike edges connecting equivalent 
vertices. All such timelike edges have the same length, and we require that — t 2 << 
min{L 2 b } in order to obtain a reasonable approximation to the time-symmetry 
condition. This completes the solution of the full two-surface initial value problem 
at a moment of time symmetry. 

4. A lattice built from prisms? 

So far the discussion of Regge calculus has been general, without the need to restrict 
attention to one particular lattice. Before we can construct a solution numerically 
we must choose a particular form of lattice to work with. In this section we discuss 
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Figure 3. A section of the axisymmetric lattice used by Dubai [ |l4[ . The relation 
to the global polar co-ordinate system is shown, together with the angles used to 
specify the extra degree of freedom. 



the application of a prism-based lattice to the axisymmetric initial value problem in 
Regge calculus at a moment of time symmetry. 

Much previous work in Regge calculus has been based on simple lattice structures, 
often built from prisms. Typical examples include the work of Lewis ]lc| , where several 
T 3 x R cosmological models were approximated with a hypercubic lattice, and Collins 
and Williams where tetrahedra were used to build each three-geometry, but the 
world-tube of each tetrahedra was not fully subdivided into four-simplices. 

The principle advantages of this prism-based approach are the ability to closely 
model the symmetries of the spacetime of interest, and the computational simplicity 
when compared with a simplicial lattice. Since there are less edges in the lattice, there 
are less hinges about which curvature is concentrated, and hence less calculations are 
required to obtain the deficit angles within the lattice. 

It is largely for these reasons that the prism approach has dominated much 
previous work in Regge calculus. The major drawback of a prism-based lattice is 
that even when all edges have been specified, the lattice is still not rigid; additional 
constraints must be provided to constrain the space of possible edge lengths. In the 
presence of high symmetry, it is possible that natural restrictions are available to 
specify the remaining degrees of freedom. This has been the case for most previous 
prism-based applications of Regge calculus. 

The simplicial approach, which constructs three-manifolds from tetrahedral and 
four-manifolds from four-simplices, is more natural if one wishes to model a complex 
spacetime. In a simplicial lattice, once all edge lengths have been specified, the 
geometry is uniquely determined. Noting that a simplicial lattice may be obtained by 
subdivision of a prism-based lattice, some of the advantages of the prism construction 
may be carried over to the simplicial lattice. 

For the particular case of axisymmetric non-rotating Brill waves, the symmetry 
suggests that a prism-based lattice may be convenient. Indeed, Dubai |Q 
approximated the pure Brill wave spacetime at a moment of time symmetry with 
a lattice constructed from the co-ordinate blocks of the continuum, a typical element 
of which is shown in figure ||. The desired axisymmetry is built into the lattice by 
demanding that there is no variation along the 0-axis. The limit as A(f> tends to zero 
is then taken to obtain a purely axisymmetric lattice. 

Dubai fixed the remaining degrees of freedom in the prisms by introducing several 
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angles on each p — z face, as shown in figure g, and specifying them in terms of the 
surrounding edges. The clever assignment chosen by Dubai allowed him to construct 
Regge solutions in reasonable agreement with the continuum solution, at least for 
low amplitude Brill waves. Dubai noted several limitations to his approach, including 
the relatively poor Regge mass estimates for "large amplitude" (a « i) Brill waves. 
Another indication of problems was the low convergence rate of the prism-based Regge 
solution to the continuum. As Dubai notes [Q , these problems arise because of the 
approximation made in specifying the angles within each prism. To obtain a better 
solution, an improved relation between the angles and prism edges is required. We 
argue that the best solution is to abandon the prisms altogether, and introduce a 
tetrahedral lattice. 

Once the axisymmetric limit {A(j) — > 0) has been taken, the major difference 
between a prism-based and a tetrahedral lattice is the absence of a diagonal edge on 
the pz faces. Each such face in a prism-based lattice is required, by construction, to 
be flat. However, the pure Brill wave metric (||) allows fluctuations across faces in the 
p — z-plane, through both the conformal factor ip(p, z) and the form function q(p, z). 
A prism-based lattice is unable to capture such variations. 

The most natural solution to this problem is to abandon the prism-based lattice 
for one constructed from simplices. In the remainder of this paper, we shall follow 
just such a path, and show that a tetrahedral lattice resolves all of the problems 
encountered by Dubai. 



5. A tetrahedral three-geometry 

In this section we construct an axisymmetric, shear-free simplicial three-geometry, and 
use it to build time symmetric initial data using Regge calculus. 

The simplicial lattice may be obtained by subdividing the basic prism shown in 
figure ^. Each block is divided into six tetrahedra, introducing three face diagonals 
and one body diagonal per vertex. The result is shown in figure f|. The body diagonal 
within the prism is denoted bijk, whilst face diagonals spanning the p and 0-axes are 
written d^t, and so forth. 

LJK 

To obtain an axisymmetric approximation in the style of the preceding section, 
we take the limit as the prism is collapsed along the </)-axis. This is a more complicated 
procedure than in the previous case ]T^ ], since we must now deal with the limiting 
form of the extra edges introduced in the simplectic approach, together with their 
associated deficit angles. 

To aid in constructing the lattice, we consider the "prototype" axisymmetric 
metric 

ds 2 = dp 2 +dz 2 + p 2 d<j) 2 . (19) 

This non-physical metric is useful to frame our discussion of the axisymmetric lattice, 
and later we will map from this "prototype" to metrics of physical interest. Figure [I] 
shows the subdivision of a co-ordinate block, such that the I, h and w edges are locally 
aligned with the p, z and </>-axes respectively. We enforce axisymmetry by explicitly 
setting 

Wijk = rij k A<t>, (20) 

and demanding that there be no variation in the edges along the 0-axis. The redundant 
j index is neglected in the remainder of this paper. 



i+lj+lk+1 




Figure 4. The rectangular prism shown in figure ^subdivided into six tetrahedra. 
This involves adding a diagonal brace to each face of the prism, together with a 
body diagonal. 



Consider the behaviour of the lattice in the limit as A(f> approaches zero. For a 
single triangle with edge lengths a, b and rAcj), in the limit as A<fi — ► we require that 
a = b. In general, the manner in which the edges a and b approach this limit is not 
uniquely defined. However, as we shall see, the assumption of axisymmetry uniquely 
determines the leading order terms in the expansions about A<f> = 0. 

Assuming that the lattice is both axisymmetric and shear- free, the p — <f> and z — <p 
faces of the block must be fiat as A(j> approaches zero, since any deviation from flatness 
introduces an asymmetry across the prism. This allows us to write an expansion in 
A<j) for the diagonal edges on these two faces. For the diagonal edge within the p — tf> 
face we have 

2 

ik 



ll+r lk r t+lk Aq> 2 , (21) 



and similarly, the diagonals which span z<j) faces of the prism may be expanded as 



{#} 



= h\ k + r ik r ik+1 A(j) 2 , (22) 



in the limit as A</> tends to zero. 

The only remaining edge which spans the y-axis is the body diagonal, bik- The 
shear-free and axisymmetric nature of the lattice also implies that the plane formed 
by the four vertices k), + 1, k), (i + k + 1) and [i + l,j + 1, k+ 1) is flat 
as A(j) tends to zero, and so the body diagonal may be expanded as 

^ = (O 2 + n fc r l+lfe+1 A0 2 , (23) 

to leading order in A(j). 

We now turn to the calculation of the deficit angles. Consider a tetrahedron 
consisting of the vertices (0, 1, 2, 3), with corresponding edges L a b, where a and b are 
vertices of the tetrahedron. The dihedral angle 6q{ 23 about the edge joining vertices 
and 1 within the tetrahedron 0123 is 

-1 

16A 12^013 



C0s6>qJ 23 — — — {L^ + 2l? m h\ z - _ -^01-^03 ( 24 ) 



"-^01^12 ^01^13 + ^02-^03 + L\ 2 L\ Z 
-^02-^13 -^03-^12/ 
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where ^.012 and A 013 are the areas of the triangular faces 012 and 013 respectively, 
given by 

16^afcc = ~-^ab _ L ac - L bc + 2L 2 ab L 2 ac + 2L 2 ab L 2 bc + 2L 2 bc L 2 ac . (25) 

Equations (|24|) and ( p5| ) are used to calculate all dihedral angles within the prism 
of figure [|. Consider, for example, the deficit angle about the edge lik- There are 
six tetrahedra hinging on this edge, so we must calculate the dihedral angle about 
lik in each of these tetrahedra, after which a power series expansion is taken about 
A(f) — 0, using equations ( pl| ) — (|23]). Retaining only the leading power in A(/>, the 
defect about Uk is 



e(hk) 



l ik ( r «fe + r i+ik - 2r i+lk+ i) + Air ik (d 2 k - h 2 +lk ) 
, l \k i r ik + n+ik - 2r lfe _i) - Air lk (d%_ t - h 2 k _^) 



where the two triangular areas constituting ap-z face are given by 



(26) 
A0, 



Afk — 7 \ + ^fk^i+ik + 2h 2 +lk d 2 k - dj k - lf k - h 2 +lk 



9/2 U2 , o/,2 J2 _ J4 _ 74 _ l2 
Zl ifc+l' t ifc ^ An ik a ik u ik l ik+l n ik' 



4 V ~ ^'+1 ^ 
with A^rjfc = r i+ i fe - r lfc . 

A similar procedure applied to the six tetrahedra about the edge r^-A0 yields 
the deficit angle 

e(r lk A(t>) = 2tt - cos — 

\ zaikiiik 

h 2 , 72 _ J2 

2^-ife/iife 

;2 , j2 _ l2 

2^zfc^i/c 

^2 , 72 _ J2 

"ifc-l + 'ifc a ik-l 

2likh-ik 1 
,1,2 ,j2 _ 72 



COS 



COS 



COS 



COS 



Zdj-ifc-i/ijfe-i 

I 2 +d 2 - h 2 

l i-lk ' "t-lfe-1 n i-\k-\ 



2ti— lfcdj-ife-i 

to leading order in A<^. The deficit angles about the vertical edge /i^ and the remaining 
diagonal edge d p ik are 

^ifc ( r ik + r ik+i - 2r i+lk +i) + A k r ik (d 2 k - l 2 ik+1 ) 



e(hik) 



ik 



h 1k ink + nk+i - 2r,_i fc ) - A k r ik (d 2 _ lk - l 2 _ lk ) 
±h lk Af_ lk 

d 2 ik ( r ik + ri+ifc+i - 2r ik +i) - Ai k r lk (l 2 k+1 - h 2 k ) 
4d ik AT k 



(27) 
A(f> 
(28) 
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d 2 lk {rjk + r l+ i k+ i - 2r l+lk ) + A ik r lk (l 2 lk - h 2 i+lk ) 
4d ik Ag 



A<t> 



where the difference operator A ik is defined by A ik r ik = r i+ i k+ i — r ik . The remaining 
deficit angles are simply 

«t) = e(4k) = <bik) = 0, (29) 



by virtue of the axisymmetric expansions, equations ( |2l| ) — (23). Note that the 
deficit about the face diagonal d p ikl the only face diagonal edge remaining after the 
limit Aip — > has been taken, is not identically zero. The fully simplicial model, 
unlike the prism method described above, can have non-zero curvature about faces in 
the pz-plane. 

Evaluating the Regge initial value equation ( |l7| ) on the simplicial lattice gives 

= hk e (hk) + h-ik e (h-ik) 

+ 2r ik A(f)e(rikA(j)) + h lk e(h ik ) + h ik -i e (h ik -i) (30) 
+ d ik e (dik) + di-ik-i e (eii-ifc-i) 

to leading order in A</> at each vertex (i, k) in the lattice. This equation may be 
evaluated once all edge lengths in the lattice are given, together with the definitions 
of the deficit angles in terms of the lattice edges given above. 

5.1. Pure Brill wave initial data 



We now consider the case of pure Brill wave initial data, as discussed in section 2.1. In 
order to solve the Regge initial value equation on the tetrahedral lattice, we perform 
a conformal decomposition, similar to the continuum approach. 

The lattice is constructed to mirror the cylindrical polar co-ordinate system in 
which the continuum Brill wave metric (||) is written. We obtain the Regge conformal 
decomposition by integrating the spacelike geodesies between vertices of the lattice, 
and assigning the geodesic lengths directly to the lattice edges. It is sufficient to take 
only the leading order terms in the expansion. 

To construct the Regge conformal decomposition for axisymmetric Brill waves, 
we must map the lattice edges from the "prototype metric" , equation ([l9]), to the Brill 
wave metric, equation (^|). Performing this mapping yields 



i Az k (31) 

2 



The notation H i+ i k indicates that the quantity H is centred between the (i,k) and 
(i + 1, k) vertices. 

The boundary conditions applied to the lattice are the same as those used to 
construct the continuum solution — equations (^) and (|l0|). They are applied using a 
second order power series expansion into the domain. On the p = boundary we take 

= Co + C lP + C 2P 2 + 0( P 3 ), 



lik 




hik = 




dik = 




r lk A(j) = 


Ipik Pi 
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Figure 5. The conformal factor ip for a Brill wave of amplitude a = 10 calculated 
using the simplicial Regge lattice. The agreement between this and the continuum 
solution in figure hi is excellent. Calculations were performed using a 601 X 601 
lattice, with the outer boundaries placed at p = 20 and z = 20. 



where the Ci are constants. Applying reflection symmetry gives C\ = 0, and evaluating 
the expansion at p = 0, Ap, 2Ap yields the condition 

0ofc = i (4^1* - ^ 2 fe) • (32) 
Similarly, the reflection boundary condition at z — takes the form 

tpio = g (4^i -ipa), (33) 

and at the inner point p = z = an expansion in both p and z into the domain 
suggests a boundary condition of the form 

000 = "001 + "010 - "011- (34) 

The outer boundary condition, equation (0) , is applied to the lattice using a centred 
finite difference approximation to the derivative. It is first converted to pure p or z 
derivatives, to give 



dip p f 1 — ip 

dp r \ r 

dip z ( 1 — %jj 

dz r 



along p = p max , (35) 



along z = z max , (36) 



and then evaluated one point in from the relevant outer boundary. 

To complete the construction, the lattice spacings Api and Azk are chosen, 
together with an initial guess for the conformal factor if> (taken to be ip = 1). The edge 
lengths on the lattice are then calculated using equations ([31]), again using equation 
(||) to specify q(p,z). Finally, the Regge-Hamilton constraint j30| ) is solved at each 
vertex, together with the boundary conditions described above. 

The Newton- Raphson method is used to linearize the system, which is then solved 
using a sparse storage biconjugate gradient algorithm. This process is repeated until 
the desired convergence criteria is reached; in our case, a fractional change in the 
conformal factor of less than 1 part in 10 -12 . Figured displays the solution to equation 
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a 


M s 


M c 


M (Alcubierre et al ) 


1 


4.7 x 10-' A 


4.7 x 10-' z 


4.8 x 10" 2 


2 


1.72 x 1CT 1 


1.72 x 1CT 1 


1.74 x 10" 1 


5 


8.77 x 1CT 1 


8.79 x 10- 1 


8.83 x 10- 1 


10 


3.22 


3.22 


3.22 


12 


4.84 


4.84 


4.85 



Table 1. Mass estimates for the different initial data sets using the Eppley 
form function q(p,z). The mass is calculated by examining the fall off of ifi 
in the asymptotic region, for each data set — the simplicial Regge (M$) and 
continuum (Mc) s olut ions. We also show the results of the previous calculations 
by Alcubierre et al ]12l , which are in excellent agreement with both the continuum 
and simplicial Rcggc solutions. 



( p0| ) with wave amplitude a = 10. The solution is in excellent qualitative agreement 
with the continuum solution shown in figure [j]. In section |6| we provide a quantitative 
analysis of the difference between the two solutions. 

We estimate the mass of the initial data shown in figures |l| and |^ by examining the 
fall off of the conformal factor in the asymptotic region. Assuming that the conformal 
factor takes the form 

M , . 

tp = 1 + - — h • ■ ■ as r — > oo, (37) 
2r 

we perform a least squares fit of the this function to each numerical solution far from 
the region in which the Brill wave is concentrated. This asymptotic form for ip is 
guaranteed by the choice of outer boundary condition — only the value of M is left 
undetermined. The mass estimates are shown in table [j] for various choices of Brill wave 
amplitude. The masses calculated for the simplicial Regge and continuum solutions 
agree remarkably well, and also show excellent agreement with previous calculations. 



5.2. Distorted black hole initial data via Regge calculus 

We now turn our attention to Regge initial data for Brill waves in a black hole 
spacetime. Since these waves may be viewed as a distortion to the Schwarzschild 
solution, it is natural in to represent them in spherical polar co-ordinates, as we found 
in section [2.2| . 

We approximate spherical symmetry in the lattice thrcc-gcomctry by aligning the 
base edges of the lattice with a spherical polar co-ordinate system. This implies, for 
example, that the z-axis in figure ^ becomes the azimuthal 0-axis. For convenience we 
also introduce an exponential radial co-ordinate rj. 

In the case of a single black hole with no Brill perturbations, the face (hik,rik) 
in figure ^ lies at constant proper distance from the throat (minimal two-surface), 
and the edges Uk stretch from one such face to the next. For the remainder of this 
section the labels i and k are taken to refer to the radial and azimuthal co-ordinates, 
respectively. This identification of the lattice edges greatly simplifies the application 
of spherically symmetric boundary conditions. 

The metric for black hole plus Brill wave initial data, equation ([ll]), is used 
to assign the base edge lengths. Integrating the geodesies between the vertices 
on each edge, and retaining only leading order terms, yields a simplicial conformal 
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(a) 



(b) 





(c) 



(d) 



Figure 6. The conformal factor for distorted black hole initial data, shown as a 
ratio of the black hole conformal factor (tp/ipbh)- The calculations were performed 
using the simplicial Regge lattice with 601 X 601 vertices, with the same choice of 
parameters as figure § (a) (1,0,1), (b) (-1,0,1), (c) (1,2,1) and (d) (-1,2,1). 



decomposition of the form 









hik = 




(38) 


dik = 


W^ lH k + .^ + A9l) 1/2 




r lk Aip = 


sin k A0, 





where we have mapped the lattice edges from the "prototype" metric, equation (|19|), 
to the distorted black hole metric. 

The inner boundary conditions are obtained as before, by combining the 
continuum inner boundary condition (|l4| ) with a power series expansion into the 
domain. The outer boundary uses a second order finite difference approximation 
to the Robin condition, equation (|l5|). The inner r\ — boundary condition is 

\ (4^i* - i>2k) , (39) 



TpOk 



and similarly, the reflection conditions on at 9 = and 6 = ir/2 are 



(40) 
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fan e = 3 (tyino-l ~ fane -2) ■ ( 41 ) 

These expansions are combined to give boundary conditions at the two "corner points" 
along 77 = 0, 

fao = "001 + fao - fax (42) 
fan e = fang-X + fang ~ fa ng -X- (43) 

These boundary conditions, together with the assignment of edge lengths given 
by equations (EjsJ), allow us to calculate and solve the initial value equation ( j3Q|) as 
before, except that we now take the pure black hole solution ([l2] ) as the initial guess. 
In the following calculations we use Bernstein's || choice of q(r],9), equation (|l6|), 
together with a choice of the free parameters (a, b, ui) and the black hole mass m. 

The solution is shown in figure |^, with the same set of parameter choices as figure 
||. We again find excellent qualitative agreement between the continuum solution and 
the simplicial Regge initial data. In the next section we will attempt to make the 
comparison more rigorous. 



6. Convergence to the continuum 



We have seen in the preceding sections that the simplicial Regge solutions show 
reasonable agreement with solutions of the continuum equations. We now make this 
comparison more quantitative. 

We have at hand two different numerical solutions for a given physical system and 
choice of the various parameters; (i) the finite-difference solution to the continuum 
equation ([l3|), and (ii) the solution of the Regge initial value equation (|30|). Since an 
exact solution is not available for comparison, the convergence analysis must be based 
solely on these equations and corresponding solutions. 

Consider the finite-difference solution fa/, to the continuum equation. It 
represents an approximation to the exact solution \P, 

fak = #<k + E lk A p + ... (44) 
where = &(T)i, Ok), A is a typical grid scale, and Eij = 0(1). We expect that p = 2, 
since a second ord er accurate approximation to equation ( pT| ) was used to generate 



fak- In section 2.2 a convergence test for the pure black hole solution confirmed that 



the continuum code is second order accurate. 

Now consider the simplicial Regge solution ipf]^, which we also expect to differ 
from the continuum solution by some amount dependent on the scale length A. By 
this reasoning we see that the Regge solution should also differ from the numerical 
solution of the continuum equations by some small amount. If we assume that 

fal = fa lk + K lk Ai + ..., (45) 

where Kij = 0(1), we can write 

^ ~ ^ k = ^*A« + (46) 
fak **fc 

where we have neglected only terms smaller than A 9 . Noting that the co-efficient of 
A 9 scales as order unity, we define the "error" between the numerical continuum and 
Regge solutions as 

fal ~ fa, 



j=0 fe=0 



4>. 



ik 



(47) 
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u -4 
10 ^ 



10 3 

100 1000 
N (number of grid points) 

Figure 7. The averaged fractional difference ejv between the simplicial and 
continuum solutions, shown as a function of the number of vertices (grid points). 
All calculations were performed on an N X N grid, for the pure Brill wave space 
with a = 10 (triangles), and distorted black hole initial data with parameters 
(1, 2, 1) (diamonds). In both cases, the fractional difference between the simplicial 
Regge and continuum solutions reduces as the second power of the grid spacing. 
From this we conclude that the simplicial Regge solution is a second order accurate 
approximation to the underlying continuum solution. 

where TV = (n n + 1) x (n$ + 1). The error ejv may be calculated directly from the two 
numerical solutions, and doing so for a variety of grid scales A yields an estimate of 
the convergence rate q. 

Figure |t] shows the behaviour of the error as the number of grid points is 
increased (hence A decreased) for both the pure Brill wave and the distorted black 
hole initial data sets. The figure suggests a value of q « 2. That is, the difference 
between the simplicial Regge solution and the numerical solution ipy, decreases as 
approximately the second power of the discretization scale. Since the numerical 
solution ipik itself differs the exact solution ^ by terms of order A 2 , we conclude that 
the Regge solution is a second order accurate approximation to 'J', the underlying 
exact solution of the continuum equations. 

7. Discussion and further work 

We have successfully used Regge calculus to construct axisymmetric initial data at a 
moment of time symmetry. Brill wave and distorted black hole initial data sets were 
obtained, and shown to agree well with the corresponding continuum solutions. 

The simplicial lattice — consisting of tetrahedra in three-dimensions — was found 
to yield solutions in excellent agreement with the continuum, and displayed none of 
the problems encountered with a prism-based Regge lattice. In particular, we observed 
that the simplicial Regge solution converged to the numerical solution of Einstein's 
equations as the second power of the lattice spacing. From this it was inferred that 
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the simplicial Regge initial data is a second order accurate approximation to the exact 
solution. 

Even for axisymmetric gravitational waves at a moment of time symmetry, we 
have found it necessary to abandon the prism-based lattice in favour of a simplicial 
approach. This indicates that the future development of Regge calculus as a 
competitive and robust alternative to standard numerical techniques will require the 
use of fully simplicial lattices. 

Work is currently underway to evolve this data set, in (2 + l)-dimensions using 
an axisymmetric adaptation of the Sorkin evolution scheme |§, , as well as a fully 
(3 + l)-dimensional evolution. 
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